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Part  I: 

Quadrature  of  HiOsc  integrals 

Suppose  that  the  kernel 

K(x,to)  £  L[[0, 1]  x  M+] 
oscillates  rapidly  for  u>  »  1,  e.g. 

•  The  Fourier  oscillator  eiuJX; 

•  The  irregular  oscillator  eiuJ9^x\  where  g  is  real; 

•  The  singular  oscillator  e'^x~y^/\x  -  y\a,  where  y  e  [0,1] 
and  a  <  1; 

•  The  Bessel  oscillator  J^(c ux),  v  e  R; 

•  The  Airy  oscillator  Ai(-cua;). 

We  are  interested  in  approximating 

I[f]  =  [  f(x)K(x,u))dx. 

JO 
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What’s  wrong  with  Gaussian  quadrature? 

Gauss-Christoffel  quadrature: 

I[f\  ~QGC[f]  =  f  4>{x)6x=  E  bif(ci)K(ci,uj), 

J0  1=1 

where  <j>  is  the  (y  -  l)-degree  polynomial  that 
interpolates  f{x)K{x,ui)  at  the  quadrature  nodes 

ci  <  c2  <  ■  •  •  <  cv  in  [0, 1]. 

In  particular:  ci,c2, . . .  ,c„  zeros  of  Pu{2x  —  1)  =>■ 
the  quadrature  is  of  maximal  order  2w.  this  is  the 
Gauss-Legendre  quadrature. 

Suppose  for  simplicity  that  K(x,  to)  =  eIUJ:  . 

For  u!  »  1  and  fixed  x  ^  0  the  value  e'UJX  is,  to  all 
intents  and  purposes,  a  random  number  on  the 
complex  unit  circle  |.~  =  1.  Therefore,  for  fixed  v 

QGC[f]  =  E  b//(q)eia;Ci  ~  0(1)  ,  a;  -  oo. 
1=1 

On  the  other  hand,  Riemann-Lebesgue  implies  that 

Jim,  Ilf]  =  0,  fe  Li  [0,1]. 
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e(i+ia;)xdx  = 


el+iw  _  1 

1  +  iu; 
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Generalising  the  Filon  method 

Instead  of  interpolating  the  integrand  f(x)K(x,u> )  at 
the  quadrature  nodes,  we  interpolate  the  values  of 
f(x)  there  by  the  polynomial  4>\ 

I[f]  ~  QF[f]  =  f  $(x)K(x,uj) 6x= 

J0  1=1 

Note  that  the  weights  depend  on  the  frequency  co. 


Filon-Legendre:  v  =  2,  c  =  [^  —  ^  + 


Just  two  quadrature  points. . . 


. . .  but  we  can  do  even  better! 

Choose  instead  Lobatto  points:  v  =  2,  c  =  [0,1]: 


The  wonder-method  is 

f(x)e'IJJXtix  ss 

where 


&i(w)/(0)  + 


eiw  _  i 


ft2(w)  = - 7- 


—  lc<; 

But  why  does  it  work  so  well? 


eICJ  —  1 

(=M*- 


Asymptotic  expansion 


Let  K(x,u> )  =  cu'-'/(-  ,  where  g  is  real  and  smooth. 
In  addition,  we  require  that  g'  7^  0  in  [0, 1] . 
Integrating  by  parts, 


1  ri  f 

/[/]  =  -/  L, 
\u  JO  o' 


1  /(x)  de'^O) 


dx 


lc 0 


gf(x)  dx 

M?(l)  /(-*-)  _  Q\ujq(0)  /(Q) 
0;(1)  ^'(0) 


-  — ^[dCZ/^O/dx] 

let; 

We  continue  by  induction.  Let 


<ro(x)  =  f(x), 

d  crm(x) 

^m+l(x)  =  - - 

dx  g  (x) 

Then,  in  the  limit, 


(X) 


/[/]-  E 

m=0 


1 

(_iu,)m+1 


<7m(0)  ig)g(Q) 
<?'(0) 


VinO)  c\u>q(l) 

g'(  1) 
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We  have 


=  /> 
ai: 

cr2: 

and  so  on. 
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// 


af  + 
s'  s' 


3#//2  - 


£ 


/4 


/- 


o"  1 

3^/'  + 


£ 


/3 


£ 


Asymptotic  quadrature 

Let 


s— 1 


Qslf]=  E 


_n  (-iw)m+i 


<7m(0)  ^iq)q(Q) 
m=0v  •  L5'(0) 

_  ^C1)  LcgCl) 

ff'C1) 

The  method  uses  s  —  1  derivatives  of  /  and 

QsU\-I[f]  ~o(w 

For  g(x)  —  x  we  have 


— s— 1 


cj  — ►  00. 


1 


<3s[/]  E  c_j  -jm+l 
m=0  ^  ' 


[/(m)(0)  -eiw/(m)(l)] 
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Filon-type  methods 

Given  nodes  ci,  <  •  •  •  <  c„  and  ni,  we 

choose  the  unique  polynomial  4>  of  degree  Zi  ni~  1 
such  that  for  all  l  — 

4>U\ci)  =  j  =  0, . . .  ,ni  -  1. 

A  Filon-type  method  is 

QF[f]  =  ![$}■ 


1.5  - 


0.5  - 


»» 


(0  CO 

u3\QF[ex]-I[ex]\ 


Blue:  Qa2\ 

Green:  Qf,  v  =  2,  m  =  n2  =  2,  ci  =  0,  C2  =  1; 

Magenta:  Qf,  v  =  3,  n±  =  713  =  2,  n2  =  1, 

ci  0 ,  C2  2  5  C3  1 • 
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THEOREM  If  ci  —  0,  cv  —  1  andn i  —  nu  —  s  then 
QF[f]  -  /[/]  ~  whenoj  — >  oo. 

Proof  Since  <2F[/]  -  [/]  =  /[<^>  -  /]  and 

<?G)(0)  =  0),  ^(l)  =  /^(l), 

for  j  =  0, 1, . . . ,  s  —  1,  the  proof  follows  from  the 
asymptotic  expansion  of /[^  - /].  □ 

Thus,  Filon  has  the  same  asymptotic  order  as  the 
asymptotic  method.  Typically  it  has  a  smaller  error 
constant,  which  can  be  further  decreased,  by  adding 
extra  nodes  in  (0,1). 

All  this  is  true  as  long  as  there  are  no  stationary  points 
of  the  oscillator  in  [0, 1],  i.e.  g'  ^  0  in  the  interval. 
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Stationary  points 

Suppose  first  that  g'(y)  —  0,  g”(y)  ^  0,  for  some 
y  g  [0, 1]  and  that  g'  ^  0  elsewhere. 


Naive  integration  by  parts  breaks  down,  since  division 
by  g'  introduces  a  polar  singularity.  An  alternative 
is  the  method  of  stationary  phase  (Cauchy,  Stokes, 
Kelvin),  except  that,  while  requiring  nasty  contour 
integration,  it  does  not  deliver  all  the  information  we 
need.  Instead,  let 


g0  (V)  =  f1  e'^Wd 
Jo 


x 


and 


/m= 


tf  N  (  N  ,  1  f1  f(x)  -  f(y)  deltJ9(x) 

j(y)M oO)  +  —  /  - jr\ - ^ - dx 

\uj  Jo  g  \  x)  tix 


—  f(y)y  oOO  + 


\(J 


jiog(l)fO)  -  f(y) 
g'O) 


_  yo)/(Q)  ~/(y) 

</(0) 


lie?  JO 


d  /(x)  -  f(y) 
tix  g'(x) 


^9(x)dx. 


11 


We  continue  by  induction.  Letting 


co  0*0  =  /(x), 

d  am(x)  -  <jm(y) 


we  have 


6x 


oo 


ff'Cx) 


m  €  N, 


0(w)  V 

ra=0  v  7 


(X) 


+  (-icu)m+1 

m=0  v  7 


J^g(O)  ^raCO) 

</(o) 

—  CTm(y) 

g'{  i) 


The  van  der  Corput  lemma  =>-  /x0(^)  =  0(w_1/  . 

Therefore,  using  the  first  s  derivatives  at  0,  y  and  1 
gives  an  asymptotic  method  with  an  asymptotic  error 

of  e>(w-s-iV 

Easy  generalisation  to  several  stationary  points  and  to 

g\y)  =  •  •  •  =  g(rHv)  =  0,  g^^Xy)  0,  r  >  l. 
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Filon  again. . . 

The  Filon  method  can  be  generalised  to  cater  for 
stationary  points.  Again,  the  idea  is  to  interpolate  to 
/  and  its  derivatives  at  [0,y\,y2,  ■  ■  ■  ,yn,  1},  where 
yi, . . . , yn  are  the  stationary  points. 
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‘Exotic’  oscillators 

Most  results  can  be  extended  to  more  ‘exotic’ 
oscillators,  for  example 

)  and  Ai(  —L OX) 

but,  clearly,  much  remains  to  be  done. 

Computation  of  special  functions 

In  the  stage  of  tentative  ideas:  using  Filon  quadrature 
for  fast  computation  of  special  functions  (e.g. 
hypergeometric  and  Bessel  functions)  for  large 
arguments,  (hopefully)  more  precise  than  using 
standard  asymptotic  formulae. 

Singular  integrals  (Hermann  Brunner,  Al  &  SPN) 
Similar  techniques  have  been  applied  to  the  kernels 

and  ela,lx_?/l  log  \x  —  y\, 

\x  —  y\a 

where  a  <  1  and  y  e  [0, 1  .  The  asymptotic 
expansion  is  more  difficult,  mainly  since  am  need  not 
be  smooth  at  y,  but  an  important  observation  is  that 
singularities  play  similar  role  to  stationary  points. 
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Multivariate  integrals 

This  is  perhaps  the  most  fascinating  chapter  of  our 
work! 


RESULT  1  Let  <T2  c  Mrf  be  a  compact  domain  with 
piecewise-linear  boundary, 


Then,  as  long  as  Vg(x)  ^  0  in  cl  Q, 


/(a Oe'^^dV 


r^j 


E 


m= 0 


1 

(-i  uj)m+d 


1 


where  each  functional  am  ,  a  periodic  function  in  u>, 
depends  just  on  /  and  its  first  m  derivatives  at  the  fcth 
vertex. 
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Consequently,  we  again  have  two  options,  both 
resulting  in  an  o(u>~s~d  quadrature:  either  truncate 
the  asymptotic  expansion  at  m  =  s  or  replace  /  by 
Hermite  interpolation  to  the  function  and  its  first  s 
directional  derivatives  at  the  vertices. 


RESULT  2  Let  Q  c  Rd  be  a  compact  domain  with 
piecewise-smooth  boundary  and  without  cusps.  Sup¬ 
pose  again  that  there  are  no  critical  points  in  cl  Q,  i.e. 
that  Vg  7^  0  there.  Then 


Q 

1 


f(x)eiLo'9(x)dV 


\uj  Jdn  ||  V^(aj)  || 


f(-X\n0n(x)TVg(x)eluJ^)as 


■T 


kWfi  ||Vff(®)||2 

where  n  is  the  outward  unit  normal 


^—rVg(x)e'^x)  dV, 


This  can  be  converted  into  an  asymptotic  expansion, 
a  Stokes-type  theorem,  “pushing”  the  integral  from  Q 
to  the  boundary.  All  this  can  be  extended  to  cater  for 
nondegenerate  critical  points  xq  e  fi,  where 

^g(x o)  =  0,  det  WTg(x0)  ^  0. 
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Part  II: 

Fredholm  equations  of  the  2nd  kind 

(Hermann  Brunner,  Al  &  SPN) 

Consider  the  problem 

£[<£](y)  =  M(y)  -  g(y),  y  e  [o,  l], 

where  g  is  given, 

JC[0](y)  =  f1  ^(x)e'^x-y^x 
Jo 

and  A  0  cr(/C). 

A  naive  approach  Cover  [0, 1]  with  the  grid 

o  =  vo  <  y\  <  •  •  •  <  vn-1  <  vn  =  i- 

Let  cj)k  w  4>(yk)  and  replace  integrals  with  Filon.  We 
obtain  a  linear  system  of  the  form 

N 

Y  =  A (j)k  -  gk,  k  =  0,  1, . . . ,  N. 

1=0 

This  will  not  work,  since  the  solution  </>  also  oscillates 
with  frequency  w,  and  this  means  that  our  asymptotics 
break  down. 
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An  alternative  We  seek  complex  numbers  Am 
and  complex-valued  functions  <frm  s.t. 


^  [0/n]  —  A  mfi 


m ■ 


Since 


K[<i>Ky)=  [Vct>(x)e'^y-X^(ix+  /Vo^e'^-^dx 

JO  Jy 


r0 

we  have 

d£[<ft](y) 

dy 


=  \UJ 


'V 


r0 


(j){x)e'^y-^  dx 


-  f1 4>(x)e'^x~y^6x 
Jy 


d2/C[<^>]  (y) _  2 


dy: 


=  (iw)  JC[(f>](y)  +  2iuxf>(y) 


But 


d  IC[4>\(y)  ,  ti2KL[4>\(y)  „ 

=  M  C y) ,  — t-2 —  =  X(P  Cy) 


dy 


dy- 
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Therefore 


A  <(>"  =  (\lo)2X4>  +  2i  u><t>. 


Let 


0(lj)  =  — 


2\uj 

“T’ 


then 


0"  +  <9Z0  =  0 


2i  _ 


Moreover, 


A^/(o)  =  d^MM  =  _ iWA<Ko), 

d  y 

=  d/C[^('1)  =  iwA^(l). 
dy 

The  condition  at  y  =  0  results  (up  to  normalization)  in 


4>(x)  —  (0  —  oj')Qi®x  +  (0  +  cv)e 

Note  however  that  0  depends  on  the  unknown 
eigenvalue  A. 


\6x 
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Using  the  boundary  condition  at  y  —  1  we  obtain 

(i 0  —  Lu)2e'®  =  {0  + 

Therefore 

(9  -  io)e^w  =  ±(6>  +  co)e~?w, 


Taking  the  plus  sign  we  obtain  the  transcendental 
equation 

...  0 

i0tan  -  =  oj, 

2 

while  the  minus  sign  yields 

•a 

10  COt-  =  oj. 

2 

The  solutions  of  these  equations  interlace:  the  first 
has  a  solution  for  Re  0  e  (2rmr,  (2 m  +  l)7r)  and  the 
second  in  Re 9  e  ((2 m  +  l)7r,  (2m  +  2)71-).  We 
observe  that 

The  real  part  of  0  behaves  like  0(m ] , 

The  imaginary  part  of  6  is  0(1)  and  small. 


20 


O' 


-O. 1 


-0.2 


-0.3 


20 


40 


60 


80 


lOO 


—0.4  - 


—0.5  - 


—0.6  - 


— O.V  - 


The  values  of  9m/( 2tt)  in  the  complex  plane  for  1  <  m  <  200. 


How  is  this  going  to  help? 


Let  (with  greater  generality) 


£[/](*/)  =  [  f(x)K(x,y)6x,  yG  [0,1]. 

Jo 

Then  the  Hilbert-Schmidt  theory  tells  us  that  K  has  a 
countable  number  of  distinct  eigenvalues  and 
eigenfunctions  {A m,4>m}-  Let 


X 


be  the  standard  real  L2  inner  product. 
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For  m  7^=  n 


A  m4>m{y)  —  [  4>m{y)K{x,y)  tix 

Jo 


Arre^JTi)  (j^n)  —  [  j  4>m(.x')4>n(.x)K(x,y)(\xdy 

t/  U  t/  U 


By  symmetry,  also 


<f>n)  =  [  f  4>m{x)4>n(<x)K{x,  y)6xdy 

JO  Jo 


and  we  just  deduce  L2  orthogonality  of  the 
eigenfunctions, 


(4>m,4>n)  =0,  n. 


Note  that  ( • ,  • }  is  not  a  positive  definite  inner  product 
it  is  complex-valued  and  it  is  entirely  possible  that 

(/,/)=  0,  f^O. 

However,  {<pm,  4>m)  ^  0,  and  that’s  all  we  need. 
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A  spectral  method 


We  expand 


OO 


f(v)  —  J2  fm<t>m(y). 
m= 0 


Therefore, 


where 


9m  — 


^  9m 

Jm  —  , 

A  Am 

m  >  0, 

0m)  

(SS  0m) 

(0m5  0m)  2(A^  —  2icj  —  u;^) 

We  thus  need  to  compute  (g,  0m)  for  a  large  number 
of  ms.  However,  if 


$m  —  C^m  i/3- 


m 


then 

(9,<M  =  (0m-u;)  f1  g(x)e^m+'am)x6x 

J  0 


+  (0m+^)  ^  g(x)e  (Prn  +  \am)xax 

Jo 
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Recall:  while  am  ~  27 rm  is  large,  \/3m\  is  small. 
Moreover,  g  is  nonoscillatory.  Therefore 

All  the  integrals  can  be  computed  very  fast  and 
accurately  by  either  the  asymptotic  method  or  a 
Filon-type  method. 


An  ongoing  challenge  is  to  generalize  all  this  to  other 
Fredholm  kernels,  e.g. 


K[f](v)=  j1  f{x)x^e'UJ\x~y\ dx 
Jo 


for  7  >  -1  and 


■  1  e\u\x-y\ 

£[f](y)  =  /.  f(x)- - — dx 


0  \x  —  y  I'r 


for  7  <E  (0, 1) 
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Pa 

So 


rt  III: 

ving  HiOsc  differential  equations 

We  commence  from  the  linear  ODE 

y'  —  A(t)y,  t>  0,  y(0)  =  y0. 

Suppose  that  its  solution  oscillates  fast,  e.g.  that  all 
the  eigenvalues  of  A  live  in  cl  C_  and  there  are  large 
eigenvalues  on  iM. 

Standard  numerical  methods  perform  very  poorly,  the 
reason  being  that  the  principal  error  term  of  a 
pth-order  classical  method  is  of  the  form 

hP+1Vp+1(y<itN)), 

where  T>pjrl  is  a  linear  combination  of  elementary 
differentials  of  order  p  +  1. 

y(t)  oscillates  with  frequency  uj  => 

\\y(p+1Xt)\\  ~  u,P+1\\ym, 

hence  ||X>p_|_1||  is  very  large! 
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An  alternative:  Change  of  variables 

To  time-step  from  tN  to  tN+1  =  tN  +  h,  set 

y{t)  =  e^~tN^Ax(t  -  tN),  t>tN, 

where  A  =  A(tN+i  .  Then 

x '  =  B(t)x,  t>  0,  x(0)  =  yN, 

where 

B(t )  =  e~tA[A(t)  -  A]etA. 

Since  e±J  oscillates  rapidly,  so  does  B{t). 

We  have  already  seen  that  high  oscillation  can  be 
turned  to  our  advantage.  The  main  idea  is  to  ‘invert’ 
the  reason  for  the  failure  of  classical  methods: 


Integrate,  don’t  differentiate! 


Specifically,  for  an  s-fold  integral  and  B{x)  a  product 
of  s  terms  from  {B(x i), . . . ,  B(xs )  , 


j  ■  ■  ■  J  B{x)dxs  ■  ■  ■  dx\ 


o 


UJ 


—s 
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The  Magnus  method  Letting 

x(t)  — 

we  have 

Q(t)  =  [  B(x) 6x 

Jo 

Jo  Jo  1  [B(x2) ,  B(x1)]dx2dx1 

1  ft  fx  1  fx  2 

+  4  /  /  /  [[B(x3),B(x2)],B(xi)]da:3dx2dxi 

ft  fX  1  fX  1 

+  T2JqL  L  [B(xz),[B{x2),B(xi)]\6xz6x2(ix1 

+  •  •  •  • 

Thus,  repeated  integration. . . . 

An  advantage  of  Magnus:  If  A(T)  lives  in  a  Lie 
algebra  g  then  y(t)  evolves  on  a  homogeneous  space 
M,  acted  upon  by  the  corresponding  Lie  group  Q. 
Using  Magnus  (with  or  without  change  of  variables) 
ensures  yN  e  M,  N  >  0. 
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An  example:  The  Airy  equation  Let  y"  +  ty  =  0. 
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Modified  Magnus,  h  =  \\yN  -  y(tN) II  ~  0(h3t~1/A) 
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Calculating  integrals:  We  use  Filon:  all  (multivariate) 
integrals  can  be  calculated  to  high  precision  using  just 
B{ 0)  and  B(h) .  As  before,  high  oscillation  helps 
computation! 

A  disadvantage  of  Magnus:  We  need  to  calculate  two 
exponentials  per  step: 

e±h,A  anc| 

This  can  be  problematic  when  the  dimension  is  large. 
However,  while  Q  is  typically  unstructured,  this  is  not 
the  case  with  A. 

Suppose  that  the  ODE  originates  in  a  semidiscretized 
PDE.  Then  often  A  is  block  Toeplitz  and  e1*1  can  be 
calculated  very  fast  by  FFT.  The  challenge  is  thus  to 
do  away  with  the  need  for  the  calculation  of  en(h\ 
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The  Neumann  method  To  avoid  the  calculation  of  the 
second  exponential,  we  abandon  Magnus  in  favour  of 
the  Neumann  expansion 


OO 

x(t)  =  ^  A fk(t)yN, 

m= 0 


where  J\fo(t)  =  I  and 


B(x i)  •  •  •  B(xm)6xm  ■  ■  ■  da,1!. 


Because  of  high  oscillation,  \\J\fm(h)\\ 
hence  very  rapid  convergence. 


Multivariate  integrals  can  be  computed  in  a  very  small 
number  of  function  evaluations,  similarly  to  Magnus 
integrals.  Again,  high  oscillation  of  B  means  that  Filon 
methods  are  very  precise. 


Numerical  results  for  Airy  are  virtually  identical  to 
Magnus,  but  the  method  comes  into  its  own  for  HiOsc 
PDEs,  e.g.  the  Schrodinger  equation. 
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Nonlinear  equations  Suppose  that 

y'  =  Mt)y  +  g(y ) 

is  highly  oscillatory.  Transforming  as  before,  but  letting 
A  =  A(tN),  we  have 

x'  —  B(t)x  +  e~tAg(tj\ r  + 1 ,  etAx). 

Let 

0'  =  t  >  0,  0(0)  =  I. 

Note  that  we  can  evaluate  0  by  either  Magnus  or 
Neumann.  Then 

x(t)  =  <P(t)yN 

+  C  -  0^Ag(tN  +  ^,e^A ®(0)df. 

Jo 

This  motivates  the  waveform  relaxation  approach, 

X[0](t)  =  yN, 

x[™+l](t)  =  <t>(t)yN 

+  f  0(f  -  Oe^g(tN  +  &  e^H(^))de 
Jo 
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Next  steps. . . 


•  Filon  without  derivatives:  Work  in  progress. 
Letting  nodes  depend  on  lo,  it  is  possible  to  obtain 
arbitrary  degree  of  error  attenuation  without  using 
derivatives; 

•  Exotic  oscillators:  For  starters,  how  to  compute 

/  f(x)  sin(wsin  nx)cix? 

Jo 

•  Multivariate  FliOsc  integrals:  What  are  the 
implications  of  the  Stokes-type  theorem? 

•  Volterra  HiOsc  equations:  The  current  approach 
doesn’t  scale  up  e.g.  to  singular  kernels; 

•  HiOsc  PDEs:  Much  further  work  required  for 
specific  PDEs,  e.g.  Schrodinger  and  Hamilton- 
Jacobi; 

•  Stochastic  DEs:  Perhaps. . . 
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